Lobster Predator Indices
Adam A. Kemberling
Updated on: 2021-10-29
Lobster Predator Indices
The idea behind the lobster predator index is to use characteristics of the fish community that prey on lobster as a biological indicator for predicting patterns in lobster abundance/biomass/recruitment.
Data Prep
The predator data used is sourced from the northeast trawl survey, conducted by the Northeast Fisheries Science Center. This survey runs every spring and fall sampling the seafloor using trawling gear towed behind a research vessel. Data from the survey includes abundance and size information for fishes found on the Northeastern US shelf.
Trawl Survey Data Prep
Data from the survey is loaded for abundance-based analyses, with some filtering done to remove stratum that are sampled less-frequently following a survey design change in 2008. Following the standard tidy-up steps, each species then has their expected length-weight biomasses estimated using published growth curve equations.
# Load clean trawl data, add length weight data
nefsc_clean <- gmRi::gmri_survdat_prep(survdat_source = "most recent")
nefsc_clean <- add_lw_info(nefsc_clean, cutoff = TRUE)Subsetting Data to the Gulf of Maine
For the predator indices we are narrowing our attention to patterns within the Gulf of Maine. This is done to relate these patterns in predators to lobster observations in the area. These include regional landings, juvenile surveys, and recruitment indices. To do this data is only kept if it is sampled from within trawl stratum that we typically associate with the Gulf of Maine
Loading Trawl Survey Strata:
# Load the strata
survey_strata <- read_sf(str_c(res_path, "Shapefiles/BottomTrawlStrata/BTS_Strata.shp")) %>%
janitor::clean_names() %>%
filter(strata >= 01010 ,
strata <= 01760,
strata != 1310,
strata != 1320,
strata != 1330,
strata != 1350,
strata != 1410,
strata != 1420,
strata != 1490)
# Key to which strata = which regions
strata_key <- list(
"Georges Bank" = as.character(13:23),
"Gulf of Maine" = as.character(24:40),
"Southern New England" = str_pad(as.character(1:12), width = 2, pad = "0", side = "left"),
"Mid-Atlantic Bight" = as.character(61:76))
# Assign Areas by Strata
survey_strata <- survey_strata %>%
mutate(
strata = str_pad(strata, width = 5, pad = "0", side = "left"),
strata_num = str_sub(strata, 3, 4),
area = case_when(
strata_num %in% strata_key$`Georges Bank` ~ "Georges Bank",
strata_num %in% strata_key$`Gulf of Maine` ~ "Gulf of Maine",
strata_num %in% strata_key$`Southern New England` ~ "Southern New England",
strata_num %in% strata_key$`Mid-Atlantic Bight` ~ "Mid-Atlantic Bight",
TRUE ~ "Outside Major Study Areas")) %>%
select(finstr_id, strata, strata_num, area, a2, str2, set, stratuma, str3, geometry)Mapping the Trawl Regions:
# Map it against the coastline
ggplot() +
geom_sf(data = new_england, size = 0.3) +
geom_sf(data = canada, size = 0.3) +
geom_sf(data = survey_strata, aes(fill = area)) +
coord_sf(xlim = c(-77, -64.5), ylim = c(34.4, 45.75), expand = FALSE) +
guides(fill = guide_legend(nrow = 2)) +
theme_bw() +
theme(legend.position = "bottom", legend.title = element_blank()) +
map_theme +
labs(subtitle = "Trawl Survey Regions")
Subset to Gulf of Maine:
nefsc_gom <- nefsc_clean %>%
filter(survey_area == "GoM")Assigning Lobster Strata
Unfortunately the trawl survey and the surveys for lobster abundances are done using a different area stratification. To perform a more direct comparison among areas we reassign the trawl stations based on where they fall within the lobster survey strata.
Loading Strata
# Lobster Strata
lobstrata <- read_sf(paste0(res_path, "Shapefiles/Statistical_Areas/Statistical_Areas_2010_withNames.shp"))
lobstrata <- rename_all(lobstrata, tolower)
# Subset the ones in GOM to plot faster
lobstrata_gom <- filter(lobstrata, id %in% c(464:467, 511:515, 521, 522, 526, 551, 561))
# Map against the Coast
ggplot() +
geom_sf(data = new_england, size = 0.3) +
geom_sf(data = canada, size = 0.3) +
geom_sf(data = lobstrata_gom, fill = "transparent", show.legend = FALSE) +
coord_sf(xlim = c(-71.5, -64.5), ylim = c(41.5, 45.75), expand = FALSE) +
map_theme +
theme(legend.position = "bottom", legend.title = element_blank()) +
labs(subtitle = "Lobster Strata in the Gulf of Maine")
Strata Re-assignment
The reassignment is done by overlaying the starting locations for the trawl stations with the survey strata. Each station is assigned based on which strata is falls within. From there we can estimated area-weighted catch rates within each of them.
Strata Re-assignment function:
# Assign statistical zone from new sf for stat zones
assign_stat_zones <- function(survdat, zone_sf, strata_col_in = "id", strata_col_out = "stat_zone", keep_NA = FALSE){
# Transfer to shorthand names
x <- as.data.frame(survdat)
out_name_sym <- sym(strata_col_out)
# use only station data for overlay/intersection
stations <- distinct(x, cruise6, stratum, station, decdeg_beglat, decdeg_beglon)
# Convert stations to sf
stations_sf <- st_as_sf(stations, coords = c("decdeg_beglon", "decdeg_beglat"), crs = 4326, remove = FALSE)
# Project to Lambert Conformal Conic
lcc <- st_crs("+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-72 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 ")
stations_sf <- st_transform(stations_sf, crs = lcc)
# Prepare statistical zones in same CRS
stratum <- st_transform(lobstrata, crs = lcc)
# rename stratum column to match desired label
names(stratum)[which(names(stratum) == strata_col_in)] <- strata_col_out
# Identify points within each polygon/strata
stations_sf <- st_join(stations_sf, stratum, join = st_within, left = TRUE)
# Don't need to convert back b/c we kept coordinates
stations_wgs <- st_drop_geometry(stations_sf)
# Keep NA's or not?
if(keep_NA == F){ stations_wgs <- filter(stations_wgs, is.na({{out_name_sym}}) == FALSE)}
# Join station assignments back into full data
out <- right_join(stations_wgs, x, by = c('cruise6', 'stratum', 'station', "decdeg_beglat", "decdeg_beglon")) %>%
mutate({{out_name_sym}} := as.character({{out_name_sym}}))
# return the table
return(out)
}Assigning Lobster Strata:
# Assign those zones!
nefsc_lobsta_zones <- assign_stat_zones(survdat = nefsc_gom,
zone_sf = select(lobstrata, id, short_name, geometry),
strata_col_in = "id",
strata_col_out = "lobster_strata",
keep_NA = FALSE)Validate by Mapping:
# make sf to check
gom_dat_sf <- nefsc_lobsta_zones %>%
distinct(decdeg_beglon, decdeg_beglat, lobster_strata) %>%
st_as_sf(coords = c("decdeg_beglon", "decdeg_beglat"), crs = 4326)
# map check
ggplot() +
geom_sf(data = gom_dat_sf, aes(color = lobster_strata)) +
geom_sf(data = lobstrata, fill = "transparent") +
geom_sf(data = new_england, size = 0.3) +
geom_sf(data = canada, size = 0.3) +
coord_sf(xlim = c(-71, -65.8), ylim = c(41, 44.5)) +
map_theme +
theme(legend.position = "bottom") +
guides(color = guide_legend("Lobster Strata", title.position = "top", title.hjust = 0.5))
Reducing to Select Strata
For this analysis we want three main strata, and one aggregate of the three. - 511
- 512
- 513
These align with our other nearshore-focused indicators.
Drop Strata with Incomplete Representation:
gom_dat <- nefsc_lobsta_zones %>%
filter(survey_area == "GoM",
lobster_strata %in% c(511:513),
#lobster_strata %not in% c(561, 522, 551, 526, 466, 467),
is.na(lobster_strata) == FALSE)Validate by Mapping:
# make sf to check
gom_dat_sf_2 <- gom_dat %>%
distinct(decdeg_beglon, decdeg_beglat, lobster_strata) %>%
st_as_sf(coords = c("decdeg_beglon", "decdeg_beglat"), crs = 4326)
# map check
ggplot() +
geom_sf(data = gom_dat_sf_2, aes(color = lobster_strata)) +
geom_sf(data = lobstrata, fill = "transparent") +
geom_sf(data = new_england, size = 0.3) +
geom_sf(data = canada, size = 0.3) +
coord_sf(xlim = c(-71, -65.8), ylim = c(41, 44.5)) +
map_theme +
theme(legend.position = "bottom") +
guides(color = guide_legend("Lobster Strata", title.position = "top", title.hjust = 0.5))
Calculate the Areas of Strata
To weight the catch rates of each stratum against one-another, they are weighted using their interior areas in square kilometers.
#### Get Stratum Area in km2 ####
# Lambert Conformal Conic
lcc <- st_crs("+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-72 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 ")
lobstrata_lcc <- st_transform(lobstrata, lcc)
# Get areas
# rename the id column to match nefsc_lobsta
# drop geometry
zone_areas <- lobstrata_lcc %>%
rename(lobster_strata = id) %>%
mutate(lobster_strata = as.character(lobster_strata),
area_m2 = st_area(lobstrata_lcc),
area_km2 = units::set_units(area_m2, "km^2"),
area_km2 = as.numeric(area_km2)) %>%
select(lobster_strata, area_km2) %>%
st_drop_geometry()
# merge in the areas
gom_weights <- left_join(gom_dat, zone_areas, by = "lobster_strata")
# Bar plot for area of each
gom_weights %>%
bind_rows(gom_weights %>% mutate(lobster_strata = "All Strata", area_km2 = sum(unique(gom_weights$area_km2)))) %>%
distinct(lobster_strata, area_km2) %>%
ggplot(aes(x = area_km2, fct_reorder(lobster_strata, as.numeric(area_km2)), fill = lobster_strata)) +
geom_col() +
scale_x_continuous(labels = comma_format()) +
labs(y = "Lobster Strata", x = "Area (km^2)", fill = "Lobster Strata")
Estimate Area-Stratified Catch
Area-stratified catch estimates are an approach that uses strata specific catch rates which are weighted by their overall size. This should be done before any species are filtered to make sure stations are not dropped by filtering out to a subset of species. At this point all the stations are in the data, and we need them to get effort correct.
Since we already reduced the potential area to these three strata, the total station count within a year is only for our area of interest. This will let us weight the three strata properly for getting the catch rates within the aggregate stratum.
Area Stratification Function:
# added some tweaks to do catch rates in km2, and expand from there
stratify_lobster_strata_catch <- function(survdat_weights, area_label, area_col){
# https://noaa-edab.github.io/survdat/articles/calc_strat_mean.html
# # Testing:
# area_label <- "lobster_strata"
# area_col <- "area_km2"
# column symbols from strings
label_sym <- sym(area_label)
area_km_sym <- sym(area_col)
#### 1. Set Constants:
# Area covered by an albatross standard tow in km2
alb_tow_km2 <- 0.0384
# catchability coefficient - ideally should change for species guilds
q <- 1
#### 2. Stratum Area & Effort Ratios
# Get Annual Stratum Effort, and Area Ratios
# The number of tows in each stratum by year
# area of a stratum relative to total area of all stratum sampled that year
# Get Total area of all strata sampled in each year (one each)
total_stratum_areas <- group_by(survdat_weights, est_year)
total_stratum_areas <- distinct(total_stratum_areas, {{label_sym}}, .keep_all = T)
total_stratum_areas <- summarise(total_stratum_areas,
tot_s_area = sum({{area_km_sym}}, na.rm = T),
.groups = "drop")
# Calculate strata area relative to total area i.e. stratio or stratum weights
survdat_weights <- left_join(survdat_weights, total_stratum_areas, by = "est_year")
survdat_weights <- mutate(survdat_weights,
st_ratio = {{area_km_sym}} / tot_s_area,
st_ratio = as.numeric(st_ratio))
# We have total areas, now we want effort within each
# Number of unique tows per stratum, within each season
yr_strat_effort <- group_by(survdat_weights, est_year, season, {{label_sym}})
yr_strat_effort <- summarise(yr_strat_effort, strat_ntows = n_distinct(id), .groups = "drop")
# Add those yearly effort counts back for later
# (area stratified abundance)
survdat_weights <- left_join(survdat_weights, yr_strat_effort,
by = c("est_year", "season", area_label))
#### 4. Derived Stratum Area Estimates
# All of the following are done based on size classes
## A. Catch per tow
survdat_weights <- survdat_weights %>%
mutate(
# Length Based :
# Catch / tow, for that year & season group
abund_tow_s = numlen_adj / strat_ntows,
# length based: biomass / tow
lwbio_tow_s = sum_weight_kg / strat_ntows,
# Not Length Based:
biom_per_lclass = (biomass_kg / n_len_class),
biom_tow_s = biom_per_lclass / strat_ntows)
## B. Catch per km
survdat_weights <- survdat_weights %>%
mutate(
# Number / km2
`abundance per km2` = abund_tow_s / alb_tow_km2,
# kg / km2
`kg per km2` = `lwbio_tow_s` / alb_tow_km2,
)
## C. Stratified Total Abundance/Biomass
survdat_weights <- survdat_weights %>%
mutate(
#### Area Extrapolations
# Number of Individuals (Expected for entire stratum, based on density)
`strat total abund` = `abundance per km2` * {{area_km_sym}},
# Sum(bodymass) of Individuals (Expected for entire stratum, based on density)
`strat total biomass (kg)` = `kg per km2` * {{area_km_sym}}
)
#### This section is for weighting stratified means, don't need here
# ## D. Area weighted catch rates
# # (catch rate, weighted by relative size of stratum)
# # stratum area : area of all stratum sampled that season
# survdat_weights <- survdat_weights %>%
# mutate(
#
# # Length Based:
# # Stratified mean abundance CPUE
# strat_mean_abund_s = abund_tow_s * st_ratio,
# # Stratified mean LW Biomass
# strat_mean_lwbio_s = lwbio_tow_s * st_ratio,
#
# # Not length based:
# # Stratified mean BIOMASS CPUE
# strat_mean_biom_s = biom_tow_s * st_ratio)
#
#
# # Extrapolated Catch Totals from area-weighted rates
# # convert from catch rate by area swept to total catch for entire stratum
# survdat_weights <- survdat_weights %>%
# mutate(
# # Length Based:
# # Total Abundance
# strat_total_abund_s = round((strat_mean_abund_s * tot_s_area / alb_tow_km2) / q),
#
# # # Two options for to estimate lw biomass | Result is the same 4/20/2021
# #
# # Option 1: Individual LW Biomass * expanded abundance at length
# strat_total_lwbio_s = (ind_weight_kg * strat_total_abund_s) / q,
#
# # # Option 2: Size specific lw biomass / tow, expanded to total area
# # strat_total_lwbio_s = (strat_mean_lwbio_s * tot_s_area / alb_tow_km2) / q
#
# # Not Length Based:
# # Total BIOMASS from the biomass of all lengths
# strat_total_biom_s = (strat_mean_biom_s * tot_s_area / alb_tow_km2) / q)
}Run Stratification:
# Run stratification
gom_strat <- stratify_lobster_strata_catch(survdat_weights = gom_weights,
area_label = "lobster_strata",
area_col = "area_km2")
# Tidy up?
# there are now two different "stratum columns" floating around
gom_strat <- gom_strat %>%
select(-c(strat_num, stratum, nafodiv, full_name)) Pull Predator Complex
The last step before estimating any of the indices is to pull out the species that prey on lobster. There are 19 species included in this list:
# List of species
lobster_predators <- c(
"atlantic halibut",
"atlantic wolffish",
"barndoor skate",
"black sea bass",
"atlantic cod",
"fourspot flounder",
"haddock",
"little skate",
"longhorn sculpin",
"ocean pout",
"red hake",
"sea raven",
"silver hake",
"smooth skate",
"spiny dogfish",
"spotted hake",
"thorny skate",
"white hake",
"winter flounder"
)# Filter
gom_predators <- gom_strat %>%
filter(comname %in% lobster_predators)
# Make another copy with both seasons included
gom_predators <- gom_predators %>%
bind_rows(gom_predators %>% mutate(season = "Both"))Summarize Predator Patterns
Depending on whether we are looking at metrics on individuals (bodysize) or across individuals (total fishes, total biomasses) there is a need to either weight/un-weight the averaging for these summaries. And an additional need for the aggregate strata to weight the contributions of each of its parts proportionally to their areas.
Individual Strata
# area of a tow
alb_tow_km2 <- 0.0384
# Independent Strata Summaries:
strata_summs <- gom_predators %>%
group_by(est_year, season, lobster_strata) %>%
summarise(
# Predator Indices (by lobster stratum):
#### Stratum Details
`stratum effort` = median(strat_ntows),
`stratum area` = median(area_km2),
#### Abundance Indices:
# Number of Individuals (Caught in Survey within stratum)
`stratum abundance` = sum(numlen_adj),
`stratified abundance` = sum(`strat total abund`),
#### Biomass Indices:
# Sum(bodymass) of Individuals (Caught in Survey within stratum)
`stratum biomass (kg)` = sum(sum_weight_kg),
`stratified biomass (kg)` = sum(`strat total biomass (kg)`),
#### Size Indices:
# Average Individual Size
# Average Individual Weight
`mean length (cm)` = weighted.mean(length_cm, numlen_adj),
`mean weight (kg)` = weighted.mean(ind_weight_kg, numlen_adj),
#### Weightings for the aggregate numbers
`aggregate area` = median(tot_s_area),
.groups = "drop") %>%
mutate(
`abundance km2` = `stratified abundance` / `stratum area`,
`biomass km2` = `stratified biomass (kg)` / `stratum area`,
)
# glimpse(strata_summs)Stratum Aggregate:
# Same idea, but the strata is all of them?
# Anywhere we took median above we need to actually pull out
agg_areas <- gom_predators %>%
mutate(lobster_strata = "Strata 511-513") %>%
distinct(est_year, season, lobster_strata, strat_ntows, area_km2)
# AND when taking averages across we want to weight them by relative area of stratum
# get stratified mean abundances/biomass
aggregate_summs <- gom_predators %>%
mutate(lobster_strata = "Strata 511-513") %>%
select(-c(strat_ntows, area_km2)) %>%
left_join(agg_areas, by = c("est_year", "season", "lobster_strata")) %>%
group_by(est_year, season, lobster_strata) %>%
summarise(
# Predator Indices (by lobster stratum):
#### Stratum Details
`stratum effort` = sum(strat_ntows),
`stratum area` = sum(area_km2),
#### Abundance Indices:
# Number of Individuals (Caught in Survey within stratum)
`stratum abundance` = sum(numlen_adj),
`stratified abundance` = sum(`strat total abund`),
#### Biomass Indices:
# Sum(bodymass) of Individuals (Caught in Survey within stratum)
`stratum biomass (kg)` = sum(sum_weight_kg),
`stratified biomass (kg)` = sum(`strat total biomass (kg)`),
#### Size Indices:
# Average Individual Size
# Average Individual Weight
`mean length (cm)` = weighted.mean(length_cm, numlen_adj),
`mean weight (kg)` = weighted.mean(ind_weight_kg, numlen_adj),
#### Weightings for the aggregate numbers
`aggregate area` = median(tot_s_area),
.groups = "drop") %>%
mutate(
`abundance km2` = `stratified abundance` / `stratum area`,
`biomass km2` = `stratified biomass (kg)` / `stratum area`,
)
# glimpse(aggregate_summs)Re-join
# Join the independent strata summaries with the aggregate strata summary
predator_summs <- bind_rows(list(strata_summs, aggregate_summs)) %>%
mutate(season = factor(season, levels = c("Spring", "Fall", "Both")),
lobster_strata = factor(lobster_strata, levels = c("513", "512", "511", "Strata 511-513")))Assess Predator Indices
To look at metrics over time we can plot each season or the aggregate of both against one another:
Plotting Function:
# Stop re-typing code
plot_pred_index <- function(pred_dat, col_select = "stratified total biomass (kg)"){
col_sym <- sym(col_select)
ggplot(pred_dat, aes(x = est_year, y = {{col_sym}}, color = season)) +
geom_line(alpha = 0.2) +
geom_point(alpha = 0.6, size = 1) +
geom_smooth(formula = y~x, method = "loess") +
scale_color_gmri() +
facet_grid(lobster_strata~season, scales = "free") +
scale_y_continuous(labels = scales::comma_format()) +
labs(x = "", y = str_replace_all(col_select, "_", " "),
subtitle = str_c("Lobster Predator Complex: \n", str_to_title(col_select))) +
theme(legend.position = "right",
axis.text.x = element_text(angle = 45, hjust = 1)) +
guides(color = guide_legend(title = "Season", title.position = "top", title.hjust = 0.5))
}Abundance
Observed
plot_pred_index(predator_summs, col_select = "stratum abundance")
Stratified
plot_pred_index(predator_summs, col_select = "stratified abundance")
Biomass
Observed
plot_pred_index(predator_summs, col_select = "stratum biomass (kg)")
Stratified
plot_pred_index(predator_summs, col_select = "stratified biomass (kg)")
Length
Observed
plot_pred_index(predator_summs, col_select = "mean length (cm)")
Stratified
#plot_pred_index(predator_summs, col_select = "stratified mean length (cm)")Weight
Observed
plot_pred_index(predator_summs, col_select = "mean weight (kg)")
Stratified
# plot_pred_index(predator_summs, col_select = "stratified mean weight (kg)")Densities
Abundance / km2
plot_pred_index(predator_summs, col_select = "abundance km2")
kg / km2
plot_pred_index(predator_summs, col_select = "biomass km2")
A work by Adam A. Kemberling
Akemberling@gmri.org